#!/usr/bin/perl -w
use strict;
use IO::File;
use Bio::SeqIO;
#use Blaster;
#use Mosaic;
use File::Temp qw/ :mktemp /;
#Copyright 2011. Jason Weirather.
#This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License v3 or any later version.
my $ivalue = 2;
my $autocolor = 'false';
my $visBlast2path = "/home/x/Code/VISBLASTWEB";
my $webPath = "/home/x/Public/visblast/work/";
# this is very functional script.  It takes a fasta file from
# standard input.  We may have an option for it to read a file.
# it will output the results of the net blast in tablular,
# dot format and or as an image file.
our %arguments;
#our $verbose;
$arguments{'eValueThreshold'} = 1;
while (scalar(@ARGV) > 0) {
  my $argument = shift @ARGV;
  if($argument=~/-(\S+)/) {
    $argument = $1;
    if($argument eq 'o') { # expect to output to a file
      if(scalar(@ARGV) < 1) { die "expected an argument for output (-h for help)\n"; }
      my $filename = shift @ARGV;
      if($filename =~ /^-/) { die "not a valied -o filename (-h for help)\n"; }
      $arguments{'outfile'} = $filename;
    } elsif($argument eq 'i') { # force input from a file
      if(scalar(@ARGV) < 1) { die "expected an argument for output (-h for help)\n"; }
      my $filename = shift @ARGV;
      if($filename =~ /^-/) { die "not a valied -i filename (-h for help)\n"; }
      $arguments{'infile'} = $filename;
    } elsif($argument eq 'e') { # force input from a file
      if(scalar(@ARGV) < 1) { die "expected an argument for output (-h for help)\n"; }
      my $evalue = shift @ARGV;
      if($evalue =~ /^-/) { die "not a valied -e evalue (-h for help)\n"; }
      $arguments{'eValueThreshold'} = $evalue;
    } elsif($argument eq 'l') { 
      if(scalar(@ARGV) < 1) { die "expected an argument for min length (-h for help)\n"; }
      my $leng = shift @ARGV;
      if($leng =~ /^-/) { die "not a valied -l length (-h for help)\n"; }
      $arguments{'length'} = $leng;
    } elsif($argument eq 'd') { 
      if(scalar(@ARGV) < 1) { die "expected an argument for -d identity (-h for help)\n"; }
      my $identity = shift @ARGV;
      if($identity =~ /^-/) { die "not a valied -i identity (-h for help)\n"; }
      $arguments{'identity'} = $identity;
    } elsif($argument eq 'g') { 
      if(scalar(@ARGV) < 1) { die "expected an argument for -g gap_width (-h for help)\n"; }
      my $gap = shift @ARGV;
      if($gap =~ /^-/) { die "not a valied -g gap_width (-h for help)\n"; }
      $arguments{'gap_width'} = $gap;
    } elsif($argument eq 'v') {
      $arguments{'verbose'} = 'true';
    } elsif($argument eq 'h') {
      print "visBlast is a perl script written by Jason Weirather 2009\n";
      print "visBlast by default reads a fasta file from stdin.  The\n";
      print "sequences in the fasta file will be blasted against one\n";
      print "another.  The header of each fasta file as usual, is\n";
      print "indicated by the '>' symbol.  The first 'word' will be\n";
      print "taken to be the ID of the sequence.  Following the ID\n";
      print "after whitespace, you may include a label for the sequence\n";
      print "if you choose.  This label can be used for using colors to\n";
      print "distinguish between nodes in the graphical output.\n";
      print "-h                  help (this screen)\n";
      print "-o filename         output text resuts to filename\n";
      print "-i filename         input fasta file from filename (default stdin)\n";
      print "-v                  verbose\n";
      print "-e evalue           choose an e value threshold for the blast (default 1)\n";
      print "-d identity         choose [0-1] a fraction of identity the sequence\n";
      print "                    is required to share (default 0)\n";
      print "-l length           chose the minimum lenght for a highscoring pair\n";
      print "                    (default 0)\n";
      print "-g gap_size         set the size of a window which is necessary to consider\n";
      print "                    two hsp's two different hits (default 30)\n";
      print "--remove_duplicates only uses one of each unique ID.\n";
      print "--dot dotfile       output the results to a dot format file\n";
      print "--png pngfile       output the results to a png format file defaults to stdout\n";
      print "                    if you are in verbose it will not print to stdout.\n";
      print "                    legend will also be printed to a file pngfile.legend.png\n";
      print "                    if you output to a file. (recommended).\n";
      print "--changeColor dotfile fastafile  \n";
      print "                                 if you want to change the color\n";
      print "                    scheme of a previously generated dot file based\n";
      print "                    on a new fasta file with all the same gene names\n";
      print "                    you can with this this command.  It skips the BLAST\n";
      print "\n";
      exit;
    } elsif($argument =~ /-(\S+)/) { # double dash word argument
      $argument = $1;
      if($argument eq 'dot') { #file name for dot file
        if(scalar(@ARGV) < 1) { die "expected an argument for --dot dot_file (-h for help)\n"; }
        my $dot = shift @ARGV;
        if($dot =~ /^-/) { die "not a valied --dot dot_file (-h for help)\n"; }
        $arguments{'dot'} = $dot;        
      }elsif($argument eq 'changeColor') { #two args
        if(scalar(@ARGV) < 2) {die "expected two arguments for --changeColor dotfile fastafile (-$argument)\n"; }
        my $dotfile = shift(@ARGV);
        my $fastafile = shift(@ARGV);
        if($dotfile=~/^-/ || $fastafile=~/^-/) { die "imporper dot or fasta filename\n"; }
        $arguments{'dot'} = $dotfile; # this is the same as the previous dot
        $arguments{'altfasta'} = $fastafile;
        $arguments{'changeColor'} = 'true';
      }elsif($argument eq 'autocolor') { 
        $autocolor = 'true';
      }elsif($argument eq 'ivalue') {
        if(scalar(@ARGV) < 1) { die "expected an argument for --dot dot_file (-h for help)\n"; }
        $ivalue = shift @ARGV;
      }elsif($argument eq 'png') { #file name for png file
        if(scalar(@ARGV) < 1) { die "expected an argument for --png png_file (-h for help)\n"; }
        my $png = shift @ARGV;
        if($png =~ /^-/) { die "not a valied --png png_file (-h for help)\n"; }
        $arguments{'png'} = $png;        
      }elsif($argument eq 'svg') { #file name for png file
        if(scalar(@ARGV) < 1) { die "expected an argument for --svg svg_file (-h for help)\n"; }
        my $svg = shift @ARGV;
        if($svg =~ /^-/) { die "not a valied --svg svg_file (-h for help)\n"; }
        $arguments{'svg'} = $svg;  
      }elsif($argument eq 'remove_duplicates') { $arguments{'remove_duplicates'} = 'true'; }
    }
  }
}


my $temp_dir = "/tmp";  # you may need to alter this directory name its where temporary files will be written
our $eValueThreshold = 1;
if(exists($arguments{'eValueThreshold'})) { $eValueThreshold = $arguments{'eValueThreshold'}; }
our $minimumHSPlength = 0;
if(exists($arguments{'length'})) { $minimumHSPlength = $arguments{'length'}; }
our $minimumHSPidentity = 0;
if(exists($arguments{'identity'})) { $minimumHSPidentity = $arguments{'identity'}; }
our $gap_window_size = 30;
if(exists($arguments{'gap_size'})) { $gap_window_size = $arguments{'gap_size'}; }
our $verbose = 0;
if(exists($arguments{'verbose'})) { $verbose = 1; }
if($verbose) {
  print "arugments: \n";
  foreach my $arg (keys(%arguments)) {
    print "  $arg\t$arguments{$arg}\n";
  }
}

my $outfile = "testoutput.txt";
our $ntAA = "nt";  # either nt or protein
# temporary file to hold the fasta file to make the database out of
#my $template = "visblasttempXXXXXXXXXXX";
#our $temp_name = mktemp($template);
my @chars=('A'..'Z','0'..'0');
my $random_string;
foreach (1..10) {
  $random_string.=$chars[rand @chars];
}
my $temp_name = "visblasttemp" . $random_string;


if($autocolor eq 'true') {
  #replace the fasta with a labeled one
  `$visBlast2path/mcl_box/execute_web_mcl.sh $arguments{'infile'} $temp_name $arguments{'eValueThreshold'} $ivalue`;
  #we need to take that temp name and use it to put the files we just generated where they belong
  #based on the infile name
  my $finalbase = "temp_1";
  my $tempbase = $webPath . $temp_name;
  my $inarg = $arguments{'infile'};
  if($inarg=~/\/([\w_-]+)\.fasta/i) { $finalbase = $webPath . $1; }
  `mv $tempbase-hmm.zip $finalbase-hmm.zip`;
  `mv $tempbase-fasta.zip $finalbase-fasta.zip`;
  `mv $tempbase-selex.zip $finalbase-selex.zip`;
}

my @lines;
my $skipBlast = 'false'; # by default we will blast
# get input either from the infile or from standard input unless we are recoloring
if(exists($arguments{'infile'})) { 
  open(INFILE, $arguments{'infile'}) or die "could not open $arguments{'infile'}\n";
  chomp(@lines = <INFILE>);
  close INFILE;
} elsif (exists($arguments{'altfasta'}) && exists($arguments{'dot'})) {
  $skipBlast = 'true';
  if($arguments{'verbose'}) { print "... we will skip blast\n"; }
} else {
  chomp(my @lines = <STDIN>);
}


our %ids;
my $uniqueName = 'unique';
my $uniqueNum = 0;
#######special mod in this loop to make non-unique names unique.
if($skipBlast eq 'false') {
  open(FASTA1, ">$temp_dir/$temp_name.fas1") or die " could not open fas1\n";
  foreach my $line (@lines) {
    ####print FASTA1 "$line\n";
    if($line =~ /^>(\S+)\s+(\S+)/) {
      my $namething = $1;
      my $class = $2; 
      if(!exists($arguments{'remove_duplicates'})) { # if this is there we will take them out later.
        if(exists($ids{$namething})) { 
          $uniqueNum++;
          $namething .= ".$uniqueName" . "-" . "$uniqueNum";
          #die "$1 duplicate entry exists. use --remove_duplicates tag to eliminate duplicate IDs\n"; 
        }
      } #otherwise we checked to see if there were duplicates.
      if(!exists($ids{$namething})) { $ids{$namething} = $class; } 
      print FASTA1 ">$namething\t$class\n";
    } elsif($line =~/^>(\S+)/) {
      print FASTA1 ">$1\n";
      $ids{$1} = '';
    } else {
      print FASTA1"$line\n";
    }
  }
  close FASTA1;
  open(FASTAIN, "$temp_dir/$temp_name.fas1") or die " could not open fas1\n";
  #foreach my $id (keys(%ids)) { print "$id\t$ids{$id}\n"; }

  my $seqin = Bio::SeqIO->new(-fh => \*FASTAIN,
                             -format => 'fasta');
  open (FASTA, ">$temp_dir/$temp_name.fasta") or die "could not open $temp_dir/$temp_name.fasta\n";

  #print "temp name is $temp_name\n";

  # store a copy of each sequence for later use
  my %sequences; #each id needs to be unique
  my @seq_array;
  my $isProtein = 0;
  my %chk_ids;
  my $sequence_count = 0;
  while(my $seq = $seqin->next_seq) {
    if(exists($chk_ids{$seq->id}) && exists($arguments{'remove_duplicates'})) { #don't print it to fasta
    } else {
      $chk_ids{$seq->id} = 1;
      $sequence_count++;
      print FASTA ">" . $seq->id . "\n";
      print FASTA $seq->seq . "\n";
      my $wseq = $seq->seq;
      if($wseq=~/[RDQEHILKMFPSWYV]+/) { $isProtein++; }
      if(exists($sequences{$seq->id})) { die "duplicate sequences\n"; }
      my $id = $seq->id;
      # strip off any part of the name following the first part to be a descriptor
      my $nsequence = $seq->seq;
      my $nseq = Bio::Seq->new( -seq => $nsequence,
                              -id => $seq->id);
      push @seq_array, $nseq;
    }
  }

  # fasta file is loaded up with the standard input now build the database
  #print "building database ... \n";
  close FASTA;
  close FASTAIN;
  my $doProtein = "F";
  #if($ntAA eq 'protein') { $doProtein = "T"; }
  if($verbose) { print "$isProtein sequences look like protein\n"; }
  if($isProtein > 1) { 
    $ntAA = 'protein';
    $doProtein = 'T'; 
  }

  #again, we can are skipping formatting the database if we are just changing color
  open(FORMATDB, "| formatdb -i $temp_dir/$temp_name.fasta -o T -p $doProtein -l $temp_dir/$temp_name.log -n $temp_dir/$temp_name");
  close FORMATDB;
  
  #print "preforming blast ...\n";
  # we have the database in a temporary directory now blast against it
  my $dbpath = "$temp_dir/$temp_name";

  #again, we can skip doing blast if we are just changing color
  if($verbose) { print "setting up blast.\n"; }
  #my $blaster = Blaster->new(dbpath => "$dbpath", type => $ntAA, eValue => "$eValueThreshold");
  if($verbose) { print "executing blast.\n"; }
  ### find out the format of the blast reports and replace this with 
  ### the non-recirpcal version of the search
  my $nums = scalar(@seq_array);
  my $prog = '';
  if($ntAA eq 'nt') { $prog = 'blastn'; }
  else { $prog = 'blastp'; } 

  chomp(my @blines = `blastall -d $dbpath -m 9 -e $eValueThreshold -p $prog -B $nums -i $temp_dir/$temp_name.fasta`);
  my @uncommented;

  #remove comments from the blast file
  for(my $i =0;$i < scalar(@blines);$i++) { 
    if(!($blines[$i]=~/^\#/)) { #not a comment
      push @uncommented, $blines[$i]; 
    }
  }
  @blines = @uncommented;
   #### there are comments in blines starting with a number sign
  #### remove all these!!!!!!!!!!

  my %values; #contains the value of the edges
  my %has_hits;
  foreach my $bline (@blines) {
    my @fields = split(/\t/,$bline);
    my @first_two = ($fields[0], $fields[1]);
    my @pair = sort {$a cmp $b} ($fields[0],$fields[1]);
    if($pair[0] ne $pair[1]) { # filter self hits
      if(!exists($values{"$pair[0]\t$pair[1]"})) {  # go with the biggest score
        $values{"$pair[0]\t$pair[1]"} = $fields[11];
      } else {
        if($values{"$pair[0]\t$pair[1]"} > $fields[11]) {
          $values{"$pair[0]\t$pair[1]"} = $fields[11];
        }
      }
      $has_hits{$pair[0]} = 1;
      $has_hits{$pair[1]} = 1;
    }
  }
  my @noHits;
  foreach my $id (keys(%ids)) {
    if(!exists($has_hits{$id})) { push @noHits, $id; }
  }

  #my %blastReports = $blaster->blast_seqs(\@seq_array,$verbose); #remove _forkable for the normal non-forked version
  # filter blast reports
  if($verbose) { print "sorting results\n"; }
  #my %filteredReports = $blaster->filterBlastReports('filter',$minimumHSPlength,$minimumHSPidentity,%blastReports);
  ### removing the option for now to identify mosaics
  #my @noHits = Mosaic->noHitIDs(%filteredReports);
  #my @multiHits = Mosaic->multiHitIDs(%filteredReports);
  #my %network = Mosaic->getNetwork(\@multiHits,\%blastReports);
  #if($verbose) { print "printing graph...\n"; }
  &printGraphFile(\%values,$outfile);
} #end skipBlast

#we can skip checking the outfile if we are changing the color
if(!exists($arguments{'changeColor'})) { # skip if we are changing color from a dot file
  if($arguments{'outfile'}) { `cp $temp_dir/$temp_name.result $arguments{'outfile'}`; 
    if($verbose) { print "saved result file\n"; }
  }
}

if(!exists($arguments{'changeColor'})) { # skip if we are changing color
  srand(time ^ $$);
  my $rand_seed = rand(100000);
  if($verbose) { print "generating dot file\n"; }
  `$visBlast2path/makeDot2.pl $temp_dir/$temp_name.result $rand_seed nolabel > $temp_dir/$temp_name.dot`;
  `$visBlast2path/makeDot2.pl $temp_dir/$temp_name.result $rand_seed label > $temp_dir/$temp_name.labeled.dot`;
} else { # if we are changing the color we need to do something here
  `./recolorDot.pl $arguments{'dot'} $arguments{'altfasta'} $temp_dir/$temp_name.result.legend > $temp_dir/$temp_name.dot`;
}

if($verbose) { # print the legend info if we are verbose
  open(LEGEND, "$temp_dir/$temp_name.result.legend");
  chomp(my @lines = <LEGEND>);
  foreach my $line (@lines) { print "$line\n"; }
  close LEGEND;
}

if(!exists($arguments{'changeColor'})) { # this won't be there if we are changing the color
  if(exists($arguments{'dot'})) { `cp $temp_dir/$temp_name.dot $arguments{'dot'}`; }
  # now that its in the temp_file.results we need to do something with it to make it a picture
  # lets convert it to dot format first.
}


if($verbose && !exists($arguments{'png'})) { 
  print "exiting because you are running in verbose mode and did not specify an output for the png file so we can't put it in stdout\n";
  cleanpath($temp_dir,$temp_name);
  die "done son\n";
}
if($verbose) { print "creating png file.\n"; }
if(!$verbose && !exists($arguments{'png'})) { # print it to the stdout
  #you won't get a legend file doing it this way.
  my @lines = `neato -Goutputorder=edgesfirst -Tpng $temp_dir/$temp_name.dot`;
  foreach my $line (@lines) { print "$line"; }
} elsif (exists($arguments{'png'})) { 
  #png to file is recommended.
  #`neato -Goutputorder=edgesfirst -Tpng $temp_dir/$temp_name.dot -o$arguments{'png'}`; 
  `neato -Goutputorder=edgesfirst -Tpng $temp_dir/$temp_name.dot -o$temp_dir/$temp_name.result.png`; 
  `neato -Goutputorder=edgesfirst -Tsvg $temp_dir/$temp_name.labeled.dot -o$temp_dir/$temp_name.result.svg`;
  `cp $temp_dir/$temp_name.result.png $arguments{'png'}`;
  `cp $temp_dir/$temp_name.result.svg $arguments{'svg'}`;
  `$visBlast2path/printLegend2.pl $temp_dir/$temp_name.result.legend`; #give it the path to dot and neato (or just dot) and then the temp file that has the legend text
  `cp $temp_dir/$temp_name.result.legend.png $arguments{'png'}.legend.png`;
  if(exists($arguments{'dot'})) { `cp $temp_dir/$temp_name.result.legend.dot $arguments{'dot'}.legend.dot`; }
  #die;
}
####die;#######
if($verbose) { print "cleaning up after png creation.\n"; }

cleanpath($temp_dir, $temp_name);
####THE END#####

sub get_nohits {
  

}

sub cleanpath {
  sleep 1;  # make sure every thing is done.  It should be already
  my $temp_dir = shift @_;
  my $temp_name = shift @_;
  if($verbose) { print "removing $temp_dir/$temp_name.*\n"; }
  #`rm $temp_dir/$temp_name.*`;
  unlink <$temp_dir/$temp_name.*>;
  return;
}

sub printGraphFile {
  my %values =%{shift @_};
  #my %network = %{shift @_};
  my $filename = shift @_;
  open(MF, ">$temp_dir/$temp_name.result");
  my %netnodes;
  foreach my $dubk (keys(%values)) {
    my @vals = split(/\t/,$dubk);
    $netnodes{$vals[0]} = 1;
    $netnodes{$vals[1]} = 1;
  }
  my @multis = keys(%netnodes);
  #print MF "\n";
  print MF "Table: visBlastOutput\n";
  print MF "Type: $ntAA\n";
  print MF "EvalueThreshold: $eValueThreshold\n";
  print MF "GapWindowSize: $gap_window_size\n";
  print MF "node1\tnode1_description\tnode2\tnode2_description\traw_score\n";
  foreach my $pair (keys(%values)) {
    my @pair = split(/\t/,$pair);
    my $node_descrip = $ids{$pair[0]};
    my $target_descrip = $ids{$pair[1]};
    my $value = $values{$pair};
    print MF "$pair[0]\t$node_descrip\t$pair[1]\t$target_descrip\t$value\n";
  }
  close (MF);
}

sub colapse_seq_inds {
  my @inds = @_; # the array of locations of conserved indicies on the query hit
  my $window = $gap_window_size; # the minimum size between conserved indicies before we call it a gap
                   # anything less than this window size will be grouped together
  my $indexOutput = "$inds[0]-";
  for(my $i = 1; $i < scalar(@inds); $i++) {
    if($inds[$i] > ($inds[$i-1] + $window)) { $indexOutput .= "$inds[$i-1] $inds[$i]-"; }
  }
  $indexOutput.= "$inds[scalar(@inds)-1]";
  #print "$indexOutput\n";
  return $indexOutput;
}

sub count_gaps {
  my $cover = shift;
  #print "$cover\n";
  my @gaps = split(/\s/,$cover);
  my $gapCount = scalar(@gaps);
  $gapCount -= 1;
  return $gapCount;
}

